#ifndef AUT_AP_2024_Spring_HW1
#define AUT_AP_2024_Spring_HW1

#include <bits/stdc++.h>
using namespace std;

namespace algebra
{
    // Matrix data structurem
    template <typename T>
    using MATRIX = std::vector<std::vector<T>>;

    // Matrix initialization types
    enum class MatrixType
    {
        Zeros,
        Ones,
        Identity,
        Random
    };

    // Function template for matrix initialization
    template <typename T>
    MATRIX<T> create_matrix(std::size_t rows, std::size_t columns, std::optional<MatrixType> type = MatrixType::Zeros,
                            std::optional<T> lowerBound = std::nullopt, std::optional<T> upperBound = std::nullopt)
    {
        MATRIX<T> matrix(rows, std::vector<T>(columns));
        if (rows == 0 || columns == 0)
            throw("“Error: Row or column is zero.” ");
        if (type == MatrixType::Zeros)
        {
            // Fill with zeros
            std::fill(matrix.begin(), matrix.end(), std::vector<T>(columns, static_cast<T>(0)));
        }
        else if (type == MatrixType::Ones)
        {
            // Fill with ones
            std::fill(matrix.begin(), matrix.end(), std::vector<T>(columns, static_cast<T>(1)));
        }
        else if (type == MatrixType::Identity)
        {
            if (rows != columns)
                throw("Error: Dimension mismatch.");
            // Fill with identity
            for (size_t i = 0; i < rows; ++i)
            {
                for (size_t j = 0; j < columns; ++j)
                {
                    matrix[i][j] = (i == j) ? static_cast<T>(1) : static_cast<T>(0);
                }
            }
        }
        if (type == MatrixType::Random)
        {
            std::random_device rd;
            std::mt19937 gen(rd());
            T lb = lowerBound.value_or(static_cast<T>(0));
            T ub = upperBound.value_or(static_cast<T>(1));
            if (lb >= ub)
                throw("Error: lowerBound is bigger than upperBound");
            else if (!lowerBound || !upperBound)
                throw("Error: No upper or lower bounds exist.");
            if constexpr (std::is_floating_point<T>::value)
            {
                // For floating point types
                std::uniform_real_distribution<T> distrub(lb, ub);
                for (size_t i = 0; i < rows; ++i)
                {
                    for (size_t j = 0; j < columns; ++j)
                    {
                        matrix[i][j] = distrub(gen);
                    }
                }
            }
            else
            {
                // For integer types
                std::uniform_int_distribution<T> distrub(lb, ub);
                for (size_t i = 0; i < rows; ++i)
                {
                    for (size_t j = 0; j < columns; ++j)
                    {
                        matrix[i][j] = distrub(gen);
                    }
                }
            }
        }

        return matrix;
    }

    template <typename T>
    void display(const MATRIX<T> &matrix)
    {
        for (size_t i = 0; i < matrix.size(); i++)
        {
            for (size_t j = 0; j < matrix[i].size(); j++)
            {
                string str = to_string(matrix[i][j]);
                cout << format("|{:^7.7s}|", str);
            }
            cout << '|' << endl;
        }
    }

    template <typename T>
    MATRIX<T> sum_sub(const MATRIX<T> &matrixA, const MATRIX<T> &matrixB,
                      std::optional<std::string> operation = "sum")
    {
        if (matrixA.empty() || matrixB.empty())
            return {};
        else if (matrixA.size() != matrixB.size() || matrixA[0].size() != matrixB[0].size())
            throw("Error: Dimension mismatch");
        MATRIX<T> matrix(matrixA);
        if (operation == "sub")
        {
            for (size_t i = 0; i < matrixA.size(); i++)
            {
                for (size_t j = 0; j < matrixA[i].size(); j++)
                {
                    matrix[i][j] -= matrixB[i][j];
                }
            }
        }
        else
        {
            for (size_t i = 0; i < matrixA.size(); i++)
            {
                for (size_t j = 0; j < matrixA[i].size(); j++)
                {
                    matrix[i][j] += matrixB[i][j];
                }
            }
        }
        return matrix;
    }

    template <typename T>
    MATRIX<T> multiply(const MATRIX<T> &matrix, const T scalar)
    {
        MATRIX<T> mat(matrix);
        for (size_t i = 0; i < mat.size(); i++)
        {
            for (size_t j = 0; j < mat[i].size(); j++)
            {
                mat[i][j] *= scalar;
            }
        }
        return mat;
    }

    template <typename T>
    MATRIX<T> multiply(const MATRIX<T> &matrixA, const MATRIX<T> &matrixB)
    {
        if (matrixA.empty() || matrixB.empty())
            throw("Error: Empty matrix exists.");
        else if (matrixA[0].size() != matrixB.size())
            throw("Error: Dimension mismatch.");
        MATRIX matrix(matrixA.size(), vector<T>(matrixB[0].size(), 0));
        for (size_t i = 0; i < matrix.size(); i++)
        {
            for (size_t j = 0; j < matrix[i].size(); j++)
            {
                for (size_t x = 0; x < matrixA[i].size(); x++)
                {
                    matrix[i][j] += matrixA[i][x] * matrixB[x][j];
                }
            }
        }
        return matrix;
    }

    template <typename T>
    MATRIX<T> hadamard_product(const MATRIX<T> &matrixA, const MATRIX<T> &matrixB)
    {
        if (matrixA.empty() || matrixB.empty())
            return {};
        if (matrixA.size() != matrixB.size() || matrixA[0].size() != matrixB[0].size())
        {
            throw("Error: Dimension mismatch.");
        }
        MATRIX matrix(matrixA.size(), vector<T>(matrixA[0].size(), 0));
        for (size_t i = 0; i < matrix.size(); i++)
        {
            for (size_t j = 0; j < matrix[i].size(); j++)
            {
                matrix[i][j] = matrixA[i][j] * matrixB[i][j];
            }
        }
        return matrix;
    }

    template <typename T>
    MATRIX<T> transpose(const MATRIX<T> &matrix)
    {
        if (matrix.empty())
        {
            return {};
        }
        MATRIX reverse_matrix(matrix[0].size(), vector<T>(matrix.size(), 0));
        for (size_t i = 0; i < matrix.size(); i++)
        {
            for (size_t j = 0; j < matrix[i].size(); j++)
            {
                reverse_matrix[j][i] = matrix[i][j];
            }
        }
        return reverse_matrix;
    }

    template <typename T>
    T trace(const MATRIX<T> &matrix)
    {
        if (matrix.empty())
            throw("Error: Empty matrix exists.");
        if (matrix.size() != matrix[0].size())
            throw("Error: Dimension mismatch.");
        int sum = 0;
        for (size_t i = 0; i < matrix.size(); i++)
        {
            sum += matrix[i][i];
        }
        return sum;
    }

    template <typename T>
    double determinant(const MATRIX<T> &matrix)
    {
        if (matrix.empty())
            throw("Error: Empty matrix exists.");
        if (matrix.size() != matrix[0].size())
            throw("Error: Dimension mismatch.");
        MATRIX<double> tmp(matrix.size(), vector<double>(matrix[0].size(), 0));
        for (size_t i = 0; i < matrix.size(); i++)
        {
            for (size_t j = 0; j < matrix[i].size(); j++)
            {
                tmp[i][j] = static_cast<double>(matrix[i][j]);
            }
        }
        if (tmp[0][0] == 0)
        {
            size_t i = 0;
            while (i < tmp.size())
            {
                if (tmp[i][0] != 0)
                {
                    swap(tmp[0], tmp[i]);
                    break;
                }
            }
            if (i == tmp.size())
                return 0;
        }
        double result = tmp[0][0];
        for (size_t i = 1; i < tmp.size(); i++)
        {
            for (size_t j = 0; j < i; j++)
            {
                if (tmp[i][j] != 0)
                {
                    double ration = tmp[i][j] / tmp[j][j];
                    tmp[j][j] = 0;
                    for (size_t k = j; k < tmp[i].size(); k++)
                    {
                        tmp[i][k] -= tmp[j][k] * ration;
                    }
                }
            }
            result *= tmp[i][i];
            if (result == 0)
                return 0;
        }
        return result;
    }

    template <typename T>
    MATRIX<double> inverse(const MATRIX<T> &matrix)
    {
        if (matrix.empty())
            throw("Error: Empty matrix exists.");
        if (matrix.size() != matrix[0].size())
            throw("Error: Dimension mismatch.");
        double value = determinant(matrix);
        if (value == 0)
            throw("Error: Singular matrix.");
        MATRIX<double> adj(matrix.size(), vector<double>(matrix[0].size(), 0));
        for (size_t i = 0; i < adj.size(); i++)
        {
            for (size_t j = 0; j < adj[i].size(); j++)
            {
                MATRIX<double> tmp(matrix.size() - 1, vector<double>(matrix[0].size() - 1, 0));
                for (size_t x = 0; x < tmp.size(); x++)
                {
                    for (size_t y = 0; y < tmp.size(); y++)
                    {
                        tmp[x][y] = matrix[x < i ? x : x + 1][y < j ? y : y + 1];
                    }
                }
                double tmp_value = determinant(tmp);
                adj[j][i] = pow(-1, i + j) * tmp_value / value;
            }
        }
        return adj;
    }
}

#endif // AUT_AP_2024_Spring_HW1